function S = rescale (s, z)

   ## usage:  S = rescale (s, z)
   ##
   ## rescale N(0,1) var. s.z to phys. units x, using PAR

   global PAR clim
   smin = PAR.scyc.fun ;

   S = s ;

   id = s.id ; z = s.(z) ;

   [nr nc] = size(z) ;
   x = nan(nr, nc) ;

   if PAR.anc
      t = date2toy(id, s.cal) ;
      x_h = [ones(rows(t),1) cos(t) sin(t) cos(2*t) sin(2*t)] ;
   endif 

   VAR = {} ; for i=1:nc, VAR = union(VAR, clim.vars(i)) ; endfor 
   for v = VAR

      v = v{:} ;
      J = strcmp(clim.vars, v) ;

      if !PAR.anc
	 eval(["x(:,J) = anom_inv(z(:,J), clim." v ".sm, clim." v ".ss) ;"]) ;
	 if all(JJ = clim.ispos(J))
	    x(:,J) = max(0, x(:,J)) ;
	 endif
	 continue
      endif

      if all(JJ = clim.ispos(J))

	 if exist([PAR.cdffit "cdf"]) == 2 && exist([PAR.cdffit "inv"]) == 2
	    for j=find(JJ)
	       jj = find(J)(j) ;
	       q = clim.(v).pp(1,j) ;
	       x0 = clim.(v).pp(2,j) ;
	       p = clim.(v).pp(3:end,j) ;
	       x(:,jj) = prb_inv(PAR.cdffit, z(:,jj), p, q, x0) ;
	    endfor
	 else
	    F = [PAR.cdffit "_inv"] ;
	    for j=find(JJ)
	       jj = find(J)(j) ;
	       q = clim.(v).pp(1,j) ;
	       x0 = clim.(v).pp(2,j) ;
	       p = clim.(v).pp(3:end,j) ;
	       x(:,jj) = feval(F, z(:,jj), p, q, x0) ;
	    endfor
	 endif 

	 ww = x(:,clim.ispos) ; if (any(I=(ww<0))(:)), ww(I) = 0 ; endif 
	 x(:,clim.ispos) = ww ;

	 x(:,J) = log(x(:,J)) ;
	 x(:,J) = x(:,J) - repmat(clim.(v).ps(1,:), rows(x), 1) ;
	 if !isempty(smin)
	    ps2 = clim.(v).ps(PAR.numhar+1:end,:) ;
	    x(:,J) = x(:,J) .* smin(ps2,t) ;
	 endif
	 x(:,J) = x(:,J) + x_h(:,1:PAR.numhar) * clim.(v).ps(1:PAR.numhar,:) ;
	 x(:,J) = exp(x(:,J)) ;

      else

	 x(:,J) = prbfit(z(:,J), clim.(v).ppinv) ;
	 if !isempty(smin)
	    ps2 = clim.(v).ps(PAR.numhar+1:end,:) ;
	    x(:,J) = x(:,J) .* smin(ps2,t) ;
	 endif
	 x(:,J) = x(:,J) + x_h * clim.(v).ps(1:PAR.numhar,:) ;

      endif 

   endfor 

   x(isinf(x)) = nan ;

   S.x = x ;
   S.names = clim.names ;
   S.ids = clim.ids ;
   S.vars = clim.vars ;
   S.lons = clim.lons ;
   S.lats = clim.lats ;

endfunction


function y = anom_inv (x, xm, xs)

   ## usage:  y = anom_inv(x, xm, xs)
   ## inverse of anom

   N = size(x) ;
   x = reshape(x, N(1), prod(N(2:end))) ;
   [nr nc] = size(x) ;

   y = x ;

   if (nargin > 2)
      y = y .* repmat(xs, nr, 1) ;
   endif 

   y = y + repmat(xm, nr, 1) ;

   y = reshape(y, N) ;

endfunction
